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Abstract 

Background: Diabetes mellitus is a group of metabolic diseases with increased blood 
glucose concentration as the main symptom. This can be caused by a relative or a total 
lack of insulin which is produced by the /i-cells in the pancreatic islets of Langerhans. 
Recent experimental results indicate the relevance of the /3-cell cycle for the 
development of diabetes mellitus. 

Methods: This paper introduces a mathematical model that connects the dynamics of 
glucose and insulin concentration with the fi-ceW cycle. The interplay of glucose, 
insulin, and p-ce\\ cycle is described with a system of ordinary differential equations. 
The model and its development will be presented as well as its mathematical analysis. 
The latter investigates the steady states of the model and their stability. 

Results: Our model shows the connection of glucose and insulin concentrations to 
the /6-cell cycle. In this way the important role of glucose as regulator of the cell cycle 
and the capability of the /i-cell mass to adapt to metabolic demands can be presented. 
Simulations of the model correspond to the qualitative behavior of the glucose-insulin 
regulatory system showed in biological experiments. 

Conclusions: This workfocusses on modeling the physiological situation of the 
glucose-insulin regulatory system with a detailed consideration of the /3-cell cycle. 
Furthermore, the presented model allows the simulation of pathological scenarios. 
Modification of different parameters results in simulation of either type 1 or type 2 
diabetes. 

Keywords: Glucose-insulin regulation, Cell cycle, Feedback loop, ODE model 



Introduction 

The term diabetes mellitus describes a group of metabolic diseases where cells, mainly 
muscle and fat cells, are not able to take up enough glucose from the blood. This can be 
due to a relative or absolute lack of insulin (cf. [1,2]). Insulin is the hormone that increases 
the permeability of the cell membrane for glucose molecules and regulates in this way the 
uptake of glucose in the cells. Therefore, a lack of insulin leads to a failure of regulation 
of glucose homeostasis and causes the main symptom of diabetes mellitus, a persisting 
increased concentration of blood sugar - in technical terms hyperglycemia. 

The more common type 2 diabetes - formerly known as adult onset diabetes - is 
characterized by insulin resistance of the target cells. Type 1 diabetes in contrast is 
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an autoimmune disease where the organism destroys the insulin producing /3-cells [3]. 
In both scenarios glucose-insulin regulation is disturbed and the adaption of /J-cells is 
insufficient to compensate for this dysfunction. 
There are three main players in the glucose-insulin regulatory system: 

1. Glucose is the energy source for the cells and is mainly obtained by carbohydrates 
in food. An elevation of blood glucose concentration is detected by the /J-cells. It 
causes them to release stored insulin molecules and to produce new insulin. 

2. Insulin is the main regulator of glucose uptake in target cells. It increases the 
permeability of the cell membrane for glucose molecules. 

3. The /J-cells are located in the islets of Langerhans in the pancreas. They store and 
produce insulin. 

The following sections describe a mathematical model for the glucose-insulin regulatory 
system that connects dynamics in the /3-cell with dynamics in the blood and the /3-cell 
cycle. The development of the model is based on the classic insulin secretion model of 
Grodsky [4], who used a packet distribution hypothesis also described by Licko [5] in 
greater detail. In these publications insulin is assumed to be stored in packets for different 
release thresholds of glucose and the main objective is to account for staircase stimula- 
tions of glucose. In our work the classic model of Grodsky is extended to variable glucose 
and adapted to the extension with insulin and glucose blood concentrations and the /J-cell 
cycle. The aim of our model is not to show biochemical or biophysical processes in detail 
but to present the core processes and interactions in a mechanistic way. The model also 
provides possibilities for extensions and consideration of additional and more detailed 
knowledge and questions. 

The paper is organized as follows. The motivation of the model, the general setup, and 
the development are presented in Section "Aim and development of the model" In Section 
"Mathematical model" a detailed description of the mathematical model is shown. The 
mathematical analysis is presented in Section "Analysis of the model" and simulations of 
the model in Section "Simulation". The results are summarized and discussed in Section 
"Discussion". 

Aim and development of the model 

Several publications [6-8] discuss the relevance of the /3-cell mass for the development of 
diabetes mellitus. Normally there is a slow turnover of /3-cells (see [9]) but the /3-cell mass 
can adapt to metabolic demands [7,8]. The concept of dynamic /3-cell mass was under dis- 
cussion for some time but is now generally accepted. Nevertheless, there is a controversy 
on the mechanisms and the precise growth factors responsible for this adaption [6,9,10]. It 
was shown that elevated glucose levels enhance /3-cell replication [11,12]. More precisely, 
Porat et al. [13] identify the glucose metabolism via glucokinase as the main positive reg- 
ulator of /i-cell proliferation. As the model in our work does not explicitly account for the 
glucose metabolism, the more general approach of glucose concentration as regulator of 
yS-cell proliferation is used. 

Mathematical models to understand the processes of the glucose-insulin regulatory sys- 
tem have a long history. Starting with the pioneering work of Bolie [14] in the 1960s one of 
the first widely used models was the minimal model developed by Bergman and cowork- 
ers (see [15,16]) in the beginning of the 1980s. Elaborate reviews of different models 
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using ODE, PDE, DDE, and integro-differential equations (IDE) are given for example in 
Makroglou et al. [17] or Boutayeb and Chetouani [18]. In the last decade, several mod- 
els dealing with the interplay of glucose, insulin, and /i-cell mass have been developed. 
For example see the work of de Winter et al. [19], Topp et al. [20], De Gaetano et al. [21], 
or the delay-model of Li et al. [22]. Other models consider particular aspects, as, e.g., 
electrical activity of /i-cells in Cha et al. [23], islet size distribution in Jo et al. [24], or glu- 
cose regulation in the whole-body system in Kang et al. [25]. Other models are designed 
to control the maintenance of normoglycemia in patients, like the compartment model 
in [26]. 

In our approach, based on the results in [6-11], instead of modeling /i-cell mass the 
whole /J-cell cycle is taken into account and plays an important role in the regulatory 
system. The main aspect of our model is the coupling of insulin storage and of insulin 
and glucose blood concentrations with the /i-cell cycle. It provides the possibility to study 
precisely the mechanism of glucose influence on the fi -cell cycle and therefore on /i-cell 
mass. The model shows the dynamics of glucose and insulin with influence of glucose 
on the /i-cell cycle. Therefore, the dynamics in the blood are directly connected with the 
mechanisms in the islets of Langerhans. 

The model analyzes the interplay of three different negative regulation feedback loops 
which live on different time scales. With elevated blood glucose concentration insulin 
release and provision is enhanced which leads to a decrease in glucose levels. Note that 
the term provision here comprises the generation of insulin, both from stored precursors 
which might dominate in the beginning and from synthesis of further insulin. 

1. The fastest feedback loop consists in a release of stored insulin immediately after 
glucose stimulus via elevated blood glucose concentrations [4]. This first insulin 
peak reaches its maximum after about three to five minutes. 

2. The second feedback loop is due to the glucose dependent enhancement of insulin 
provision. This has a visible effect after about 10 minutes [4]. 

3. The slowest feedback loop consists of the enhancement of the /i-cell cycle by 
glucose. If the first two reactions of the system are not sufficient to end 
hyperglycemia, the blood glucose concentration remains at an elevated level. This 
mild hyperglycemia results in an enhancement of the /i-cell cycle leading to more 
/i-cells which in turn can produce further insulin (see [9]). 

There are different processes that increase /i-cell mass via cell number [8,9]. Besides repli- 
cation of existing cells, there is also neogenesis by transdifferentiation and stem cells. As 
an assumption in our paper, based on publications [27-29], the adaption of /i-cell mass is 
managed by replication only. 

Figure 1 shows a schematic concept of the model. For simulation of the regulatory 
system the model is stimulated via elevated blood glucose level. 

Besides the uptake of glucose through food, there is also a production of glucose by 
the organism itself incorporated into the model with a constant production rate. The 
glucose stimulus induces the organism to release stored insulin and to enhance insulin 
provision. Insulin is stored in packets for different release thresholds [4]. The storage 
is filled through glucose dependent insulin provision and cleared at a constant secre- 
tion rate. Secreted insulin regulates the uptake of glucose in muscle and fat cells. Besides 
the immediate release and enhanced provision per cell, the model also accounts for the 
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Figure 1 Schematic concept of the model. Glucose is given to the system at a constant production rate 
mainly by the liver. Elevated blood glucose levels lead to immediate release of stored insulin and an 
enhanced insulin provision. Also, glucose influences the transition rate between phases Gi and S of the cell 
cycle. Insulin regulates the uptake of glucose in target cells. The molecules are stored in packets with 
different release thresholds. These packets can be redistributed within the storage. 



slowest regulation feedback loop, i.e. glucose influencing the /i-cell cycle. The replication 
of /3-cells eventually leads to the provision of more insulin. 

The aim of our work is to describe the three different feedback loops in one model and 
to provide a basis for understanding and explanation of the mechanisms in the glucose- 
insulin regulatory system. 

The different parts of the model will be described in detail in the following section. 

Mathematical model 

/J-cell cycle 

The mathematical model consists of three parts where the first one is the /J-cell cycle. 
The model accounts for three phases of the cell cycle [30]: 

1. The Gi-phase is a growth phase where the cell prepares for synthesis. A basic 
assumption of the model is that the functioning /3-cell mass lies in this phase [31]. 

2. The S-phase is the synthesis phase where DNA replicates. 

3. The G2/M-phase is the premitosis and mitosis phase where the nuclear division 
takes place. 

Biological experiments concerning the cell cycle are often done by flow cytometry. 
This method measures the DNA content in the different phases and can not distin- 
guish between phases Gi and M that have the same DNA content. Therefore, both are 
combined to one G2/M-phase. 
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Figure 2 /i-cell cycle. Cell cycle of the (8-cells with three phases (Gi,S, GjjM) and transition rates (p1.p2.p3). 
The apoptosis rate is given with parameter pt,. The term pi (1 + psG(f)) describes a linear influence of 
glucose on the transition rate pi with influence factor ps. 



As it is shown in Figure 2 glucose modifies the transition rate from G\- to S-phase. 
Several authors claim this transition to be an important checkpoint for the regulation of 
the /i-cell cycle [32]. The linearity of the glucose influence, p\[ 1 + psG], is a simplifying 
assumption in our work. 

In the model, the transition rates pi,p2,P3: the apoptosis rate pn, and the influence 
factor of glucose ps are considered. All parameters can be found in Table 1. 

An important contributing factor to the /i-cell dynamics is glucose toxicity described by 
Unger et al. [33]. This term refers to the wide range of harmful effects of chronic hyper- 
glycemia leading to chronic oxidative stress after the onset of diabetes, including damages 
to the pancreatic islet /i-cell (cf. [34,35]). As a consequence of hyperglycemia lipid toxicity 
may additionally damage /S-cells. This effect is called glucolipotoxicity and is described in 
[36]. To account for the effects of glucose toxicity and glucolipotoxicity on insulin secre- 
tion a glucose dependent apoptosis rate can be incorporated to the cell cycle. In the actual 
version of the model the role of glucose toxicity is omitted for the sake of simplicity and 
the apoptosis rate is constant. As an assumption in the /i-cell cycle model, apoptosis and 



Table 1 Model parameters 





Value 


Definition 


Pi 


6.0594 x 1 0~ 5 i 


transition rate 61 S 


Pi 


4.9861 x lO^i 


transition rate 5 -> G2/M 


P3 


8.9444 x lO -4 ^ 


transition rate G2/M -*■ 61 


Pa 


3.3194 x10- 4 ^ 


apoptosis rate 


Ps 


0056 


influence factor of glucose 


P6 


0 3 mg 
U D 100 ml min 


rate of glucose production 


P7 


0 003 m 


glucose effectiveness 
at zero insulin 


Ps 


360x10-3^ 


insulin sensitivity 


P9 




secretion rate of insulin 


PlO 




decay rate of insulin 


Pl1 


0 0337 m 


rate of provision increase 


Pl2 


1.72 x 10~ 9 mg 


average insulin amount per /3-cell 


bv 


3.33 ml 


blood volume of 35g mouse 


f 


0.5 


proportionality factor 


h 


5 


Hill coefficient 


Po 


186 - 506 T0TO 


Michaelis constant 


Xmax 


1.65 x 10~ 3 mg 


maximal amount of insulin 
per pancreas 


k 


3.3 


Hill coefficient 


C 


1 AQ 19. m 9 
>^-'° 100ml 


Michaelis constant 



Parameters of model (9) used for the simulations shown in Figures 5 and 6 in Section 'Simulation'. The parameter values are based 
on biological experiments of [4,20,31,32]. 
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the only under some pathological conditions relevant necrosis are subsumed in the rate 
P4 called apoptosis rate [37]. 

The /i-cell cycle is modeled as a three compartment model as it is common in cell cycle 
modeling (cf. [38]): 



The constant 2 in the first equation accounts for cell division in the transition from G 2 /Al- 
to Gi-phase. In the physiological case for adults the /i-cell cycle is very slow (see [9]) 
but has the capability of dynamic adaption to metabolic demands [7]. In model (1) glu- 
cose influences the transition rate p\ from G\- to S-phase and is able to regulate the 
/i-cell cycle. This is the case if glucose triggers the system and neither the immediate 
release of stored insulin nor the enhanced insulin provision is able to lower blood glucose 
concentration. Then a high level of glucose forces the cell cycle to accelerate. 

Glucose and insulin concentration in the blood 

The dynamics of blood glucose concentration are based on the model of Topp et al. [20]. 
There, the change in blood glucose concentration is modeled as the difference between 
production and uptake of glucose, 

G(t) = production - uptake = p 6 - \pi + PsHt)] G(t). (2) 

The net rate of glucose production is represented by a constant production rate p 6 . This 
rate is the difference of an intrinsic glucose production, mainly by the liver, and glucose 
concentration independent uptake of glucose. The latter consists mainly of glucose uptake 
by the brain and other nervous tissues which is assumed to be constant in our model. 
The uptake of glucose in other tissues consists of two processes dependent on the glucose 
blood concentration. One is an insulin independent uptake represented by the parame- 
ter of glucose effectiveness pj. The other process is an insulin dependent glucose uptake 
mainly by muscle and fat cells which is influenced by insulin sensitivity ps and depends 
on blood insulin concentration/. The parameters are listed in Table 1. 

Similarly, the dynamics of blood insulin concentration are modeled as secretion minus 
degradation, 



where secretion consists of the secreted amount of insulin molecules from the /J-cells, 
P9X1, that has to be considered with respect to the blood volume bv of the organism. 
Variable X\ will be discussed in detail in the following section. Degradation of insulin is 
modeled with a constant decay rate p\§. 

Insulin storage 

In 1972, Grodsky published a packet distribution hypothesis for insulin granules in an 
insulin storage [4]. In this work the storage is modeled with no dynamic connection to 
the remaining glucose-insulin regulatory system. For stimulation of the system, different 
glucose functions were considered, for example single- or two-step constant stimulation, 
staircase stimulation, or ramp functions of glucose concentration. 



Gi(t) = 2p 3 G 2 /M(t) - [ Pl [1 +p s G(t)] +P4] G x (t), 
S(t) = Pl [1 + P5 G(t)] Gi(t) - p 2 S(t), 
G 2 /M(t) = p 2 S(t) - p 3 G 2 /M(t). 



(1) 



1 

Kt) = t P9Xi(f) - p 10 I(t), 



(3) 



Gallenberger etal. Theoretical Biology and Medical Modelling 2012, 9:46 
http://www.tbiomed.eom/content/9/1/46 



Page 7 of 22 



In our approach the model of Grodsky is incorporated into a model of the glucose- 
insulin regulatory system including the adaption of /3-cell mass to metabolic demands. 
Although the publication of this insulin secretion model is several years ago, the packet 
distribution hypothesis still finds application, for example in the work of Overgaard 
et al. [39]. There, it is included into a mathematical model for insulin secretion applied 
to IVGTT and OGTT data. Also Pedersen et al. presented an updated version of the 
hypothesis for oral minimal models of insulin secretion in [40] and Tsaneva-Atanasova 
described insulin secretion in a general context of mechanisms of cell secretion [41]. 

As the knowledge about /J-cell biology has increased since 1972 some reinterpreta- 
tion of the assumptions in [4] are appropriate. In the original work there is no clear 
definition of what the packets are. They could be interpreted as insulin containing 
granules within a cell with different sensitivities for glucose induced release. Although 
granules clearly are present and variability of sensors for signals is a common feature 
in biology, there is no direct experimental proof, yet, to our knowledge. An alternative 
interpretation would be variability of sensitivity on cell level. In fact, Jonkers and Hen- 
quin [42] show a sigmoidal distribution of active /i-cells. A combination of both aspects 
(and possibly more unknown factors) may be the most probable explanation for the 
experimentally found dose-response curves [43]. Therefore, in the following the pack- 
ets are interpreted as pancreatic /i-cells that can be active or inactive. Several recent 
models for insulin secretion considering Ca 2+ -evoked exocytosis, the cAMP amplify- 
ing pathway, and actual results on granule dynamics are available (see e.g. [44-46]) but 
for the qualitative conclusions in our work Grodsky 's basic model of insulin secretion 
is sufficient. 

Our first modification of the insulin secretion model [4] is the adaption to time depen- 
dent glucose concentrations. In doing so, it obtains a wider field of application and gets 
connectable to the dynamic regulatory system. In a second step the insulin storage is 
incorporated with blood glucose and insulin concentrations as well as with the /i-cell cycle 
to form a complete model. It is assumed that insulin is stored in homogeneous packets in 
the storage which is shown in Figure 3. There are different release thresholds 8 (in units of 
glucose concentration) for the packets and the major part of them is stored at lower values 
of glucose concentration, between 100 and 200 j^fg. These lower values are often reached 
in an organism so that there is the need for sufficient insulin to react to these impulses. 

The initial packet distribution f(0,O) f° r threshold 0 at time t = 0 is crucial for the 
storage model. In every time step the different processes within the storage try to achieve 
this initial distribution at least qualitatively. 

The storage is filled through a glucose dependent insulin provision factor P, and insulin 
is released with constant secretion rate pg. Besides that, there is a redistribution process 
where packets that were not needed so far are redistributed to qualitatively reestablish 
the initial packet distribution $j(0, 0). 

Using the packet distribution the storage can be divided into two compartments. The 
releasable amount of insulin X\ contains the packets with threshold value 9 below the 
actual blood glucose concentration G, i.e., 



G(t) 




o 
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e 




100 ml 



Figure 3 Insulin storage and packet distribution. The storage is filled through glucose dependent insulin 
provision P(t) and cleared at a constant rate pg. Insulin is stored in packets at different release thresholds 9. A 
special characteristic of the storage is the redistribution process of packets shown by the vertical arrows. The 
figure is adapted from [4]. 



These packets can be released from the /J-cell for glucose concentration G. The second 
compartment in the storage pool is the non-releasable amount of insulin X 2 of the packets 
with threshold value 6 above the actual blood glucose concentration G, i.e., 

00 

X 2 it)= j %(6,t)dQ. 

G(t) 

These packets can not be released for glucose concentration G but are involved in the 
redistribution process. 

As, in our work, the insulin storage is connected to the /J-cell cycle, changes in the 
insulin producing /J-cell mass influence also the initial distribution £ (6, 0) as it is depen- 
dent on the maximum amount of insulin per pancreas, X max . The insulin producing /J-cell 
mass can change in every time step t and thus X max is time dependent, too. Therefore, 
this distribution is called target distribution and is denoted by f *(6*, t). Consequences of 

this modification are shown later in this section. 

The insulin storage is modeled as a three compartment model with X\, X 2 , and an 
equation for the dynamics of the provision factor P. A detailed version shows the different 
processes within the storage according to the target distribution: 

X\(t) 

Xiit) = t(G(t),t)G(t) -psXtf) 

G(£) oo G(£) 

+f J f(9,t)de[Xi(t)+X 2 (J:)] -f j $*(0,t)deXi(f) +f j t(0,t)d6P(t), 
o oo 

^ (t)= ^W^ r(G(0 '° 6(t) ^ 

OO oo oo 

+/ J ?(o,f)do[x 1 <f)+x 2 (t)] -f j i-*(e,t)dex 2 {t) +f j $*(.o,t)d$p(f), 

G(f) 0 G(£) 

P(t)=pn[Poo(t)-P(t)]. 



For details concerning the derivation of equations (4) see the Appendix "Derivation of glu- 
cose dependent insulin storage dynamics" section. The first terms in the equations for X\ 
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and X 2 correspond to the time dependence of the integral limits (chain rule) which is an 
expansion of the original insulin secretion model. They describe the influence of changing 
glucose concentration on the separation of the insulin storage into compartments X\ and 
X 2 (see Figure 4). This influence depends on the actual glucose concentration value and 
the change in glucose concentration with time t. 

System (4) shows the production and redistribution process in detail according to the 
target distribution %*(6, t). Their contribution to the dynamics of the single compart- 
ments are modeled with proportionality factors / and /, respectively. 

To simplify the representation the integral functions are expressed in terms of general 
transition functions. The simplified version of the model is given in the following: 



Xi(t) = 



Z*(G(t),t)G(t) 



~r i u 2 (t) 



U\{t) -P9 



Xi(t) + u 2 (t)X 2 (t) + Ui (t)G x (t)P(t), 



X 2 <f) = ui(t)X 1 (t) 



S*(G(t),t)G(t) 



f~ l ui{t) 

P(t) =/>!! [Poo (*)-*(*)]. 



U2(f) 



X 2 (t)+U4,(t)Gi(t)P(t), 



(5) 



The parameters of the model are given in Table 1. The glucose dependent transition func- 
tions U\, . . . , K4 describe provision and redistribution processes according to the target 
distribution. They are given in greater detail in Appendix "Transition functions" section. 

The function P^o models the glucose dependent steady state of insulin provision. To 
express the delay in the effect of insulin provision, Poo is modeled as a Hill function 



Poo(t) = 



G(t) h 



ph 
1 0 



G(tr 



with h being the Hill coefficient and Pq the Michaelis constant. 




0 50 - 100 150 200 250 300 350 400 450 500 



glucose G and threshold 8 in mg/100 ml 

Figure 4 Initial distribution in the storage and compartments Xi andX2. Secreted insulin X(G,0) 
(dashed line) is the whole amount of insulin releasable by the pancreas at glucose concentration G. It is the 
integral over the initial packet distribution ?(0,O) (solid line). Xi is the releasable and X2 the non-releasable 
amount of insulin. In the plot they can be identified as areas under the solid curve 0) ranging from 
8 =[0, G] and from 9 =[G, 00), respectively. The figure is adapted from [4]. 
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An example for a possible initial distribution in case of constant glucose is given in 
[4]. For constant glucose concentration G experimental results [47,48] show the following 
course of insulin release X(G,0) in the early phase (dashed line in Figure 4). This evolution 
is described with a sigmoid function 

G k 

X(G,0) -Xmax ck + Gk , 

where X max is the maximum amount of insulin per pancreas, k the Hill coefficient, and 
C the Michaelis constant. With X(G,0) being the amount of released insulin at constant 
glucose concentration G, the initial distribution f(0,O) is expressed as the derivative of 
X(G, 0). To be precise, 



X(G,0) , = j 2U>A))<1» (6) 



' C k + G k 



d - kC k 0 k ~ 1 

H0,0) = -X(0,0)=X max ^ k - ¥V2 . (7) 



These expressions change in the case of time dependent glucose and in connection with 
the yS-cell cycle that results in variable /3-cell mass. First, the maximum amount of insulin 
per pancreas X max is not constant anymore but depends on the /J-cell mass G\ at time t, 

X max (t) = puGi(i). (8) 

The parameter pi2 is interpreted as the amount of insulin per /3-cell. One of several pos- 
sibilities to determine this parameter is the quotient of the maximum amount of insulin 
per pancreas and the initial amount of /J-cells, i.e., 



Pl2 



Gi(0)' 



where Gi(0) serves as a normalization. Then equation (8) determines for every amount 
of /J-cells G\ at time t > 0 the corresponding maximum amount of insulin. Therefore, 
expressions (7) and (6) now read as 

G(t) 

Gitr C 
X*{G{t),t)=X max {t) —, = j H*{O,t)d0 



C k + G(t) k 

o 

d kC k 0 k ~ 



l 



ne,t) = -x*(e, t )=x max ( t)(ck + ek)2 , 



where f * (6, t) is now called target distribution. 
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Complete model 

In the previous subsections the partial models have been developed and explained in 
detail. Based on this outline the complete model can now be formulated: 



Gi(t) 
S(t) 
G 2 /M(t) 
G(t) 

/(*) 
*i(t) 



: 2p 3 G 2 /M(t) - [pi [l+/» 5 G(t)] +/> 4 ] Gi(t), 
= /»i [l+/> 5 G(t)] Gi (t)-p 2 S(t), 

■ PiS(t) - p 3 G 2 /M(t), 

■ P6 - [P7 + p&I(t)] G(t), 



(9) 



/*9*i(t)-/>io/(t). 

-u 1 (t)-p 9 



bv 

t(G(t),t)G(t) 



f- l u 2 (t) 



Xi{t) + u 2 (t)X 2 (t) + u 3 (f)Gi(t)P(t), 



x 2 (t) = Ul (tyxi(t) 



t(G(t),t)G(t) 



Pit) 



f~ l u 2 {t) 
pn [Pooit) - P(t)] . 



+ u 2 (t) 



X 2 (t) + in(t)Gi(t)P(t), 



The first three equations describe the dynamics of the /3-cell cycle with its three phases. 
The fourth and the fifth equation show the dynamics of blood glucose and insulin con- 
centration, and the last three equations present the insulin secretion model. The partial 
models are connected through several entities. 

1. Insulin I influences the glucose dynamics via insulin dependent uptake in target 
cells. Secretion of insulin consists of the releasable amount of insulin molecules 
pgXi in relation to the blood volume bv. 

2. Glucose G plays an important role in regulating the processes within the insulin 
secretion model. It regulates the provision of insulin and defines the compartments 
of the insulin storage via the target distribution. It also contributes to the 
redistribution process. Furthermore, glucose regulates the /J-cell cycle via the 
glucose dependent transition rate from G\- to S-phase. 

3. The /S-cell mass G\ determines the capacity of insulin provision. 

In Section "Simulation" the behavior of solutions of the complete model can be seen. 
There, the model is simulated in the physiological case as well as in an experimental 
situation. 



Analysis of the model 

In this section a basic mathematical analysis is presented to achieve a better understand- 
ing of the model behavior. First, positivity of the solution is shown. Then the analysis 
focusses on steady states and their stability to explain the asymptotic development of the 
solution. 



Positivity of solutions 

Our model of the glucose-insulin regulatory system describes the dynamics of biological 
quantities, e.g., cell numbers, concentrations, and mass. Naturally, these quantities are 
positive and therefore positivity of the solution is a desired characteristic of the system. 
A necessary and sufficient condition for the existence of positive solutions is given in the 
following corollary. 
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Corollary 1. ([49]): For every t 0 > 0 and initial value x 0 = (x®, . . -,x° m ) € = {x e 
~R m : Xi > 0, i = 1, m}, each component xu i = 1, . . . , m, of the solutions to 

Xi(t) = fi(t,X\, . . . ,Xwi), Xi(to) = X^ 

is positive if and only if 

fi(t,xi, . . 0,*j+i, . . . ,x m ) > 0, (10) 
i = 1, . . . , m, for all t > 0 and x e M+. 

Using the condition from Corollary 1, positivity of the solution of system (9) can be ana- 
lyzed. The parameters of the model and the initial values are non-negative (see Tables 1 
and 2, respectively). Therefore, it can be shown that condition (10) holds for all positive 
values of the variables. Summarizing, positivity of the solution of the presented model for 
the glucose-insulin regulatory system is guaranteed. 

Steady states 

As it can be seen in Section 'Simulation', the graph of the solution suggests a steady state 
behavior of the system. The investigation of the steady states is based on the /J-cell cycle 
model (1). 

We determined the only two steady states for the /S-cell cycle. There is a trivial steady 
state where no cells are present, 

G\ = S* = G 2 /M* = 0, (11) 

and there is another steady state where the number of apoptotic cells equals the number of 
new cells. This occurs if the apoptosis rate equals the transition rate from G\- to 5-phase: 

P4 = Pl (l+p 5 G) O G = ^^i. (12) 

PlP5 

Equation (12) results in a fixed value for glucose, G, that can be modified via the influence 
factor p5. G is a crucial threshold for the development of the cell cycle as it determines 
the values of glucose concentration leading to /J-cell mass increase or decrease. 

The apoptosis rate p4 is greater than the transition rate from G\- to 5-phase if the actual 
value of glucose concentration G t is lower than G: 

p 4 > pi(l+p 5 G t ) O G>G t . (13) 

In this case there are more /i-cells dying than dividing per time step and in consequence 
the total amount of /J-cells in Gi-phase is decreasing. 



Table 2 Initial values 



Variable 


Initial Value 


Definition 


6,(0) 


958000 


number of cells in 61 -phase 


S(0) 


14000 


number of cells inS-phase 


6 2 /M(0) 


28000 


number of cells in G2/M-phase 


6(0) 


200 mg 
zuu 100 ml 


blood glucose concentration 


1(0) 


001 mg 
u ' ul 100 ml 


blood insulin concentration 


X,(0) 


0.0012 mg 


releasable amount of insulin 


X 2 (0) 


0.0005 mg 


non-releasable amount of insulin 


P(0) 


0 


provision of insulin 


Initial values and definition of variables in model (9) used for the simulations shown in Figures 5 and 6. 
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In contrast, the apoptosis rate p^ is lower than the transition rate from G\- to S-phase 
if the actual value of glucose concentration Gt is greater than G: 

< Pl (l+p 5 G t ) & G < G t . (14) 

Fewer /S-cells are dying than dividing per time step and in consequence the /J-cell mass in 
Gi-phase is increasing. 

In summary, the further analysis of the steady state behavior results in only two fixed 
points for positive values of the variables. The two steady states (11) and (12) of the cell 
cycle model determine the two steady states of the whole model. The first one is the trivial 
steady state 

F* = (o,0,oA 0,0,0, Poo (15) 

with all cell numbers in the three phases equal to zero. As it can be seen in Fj", glucose 
concentration is never equal to zero due to constant production of glucose by the liver. 
Thus, the steady state of the provision factor, Poo, is also not equal to zero but a fixed value 
depending on G* = j^- 

The second steady state is driven by the threshold value G determining a non-trivial 
steady state 

F| = (G* V S* ,G 2 /M* ,G* J* ,X* lt Xl,P*) . (16) 
The values for the different variables can be given explicitly: 

Gt = S* = G 2 /M* = P -±S* - PiP9 ^ 



PnfP* PiPxifP* ' P3 ViPxifP* 

G* = G = m ~ Pl i* = P6 ~ P?G * 



PiPs PsG* 

xr = b -^, xi = x l [ P9 [l + G*-^]-f], p*= G * h 



The steady state Ff is reached after a glucose stimulus to the regulatory system. The 
variables tend to these values if no further impulse or modification to the system is 
following. 

Stability 

The stability of the two steady states Fj" and can be investigated by computing the 
Jacobian matrix of these fixed points. This analysis is based on the system parameters in 
Table 1. 

In summary, it can be shown that there are two types of steady states for the system. 

1. The trivial steady state with cell numbers equal to zero, 

F* = (o, 0, OA 0,0,0, Poo 

V Pl \P7 

is an unstable fixed point. 

2. The non-trivial steady state, 

Ff = (G* V S*,G 2 /M*,G*,I*,X$,X 2 *,P*), 

is a stable fixed point reached by the system after some time without any influence 
from the outside. 
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The behavior of the model according to the mathematical analysis will be illustrated in 
the following section using simulations of the complete model. 

Simulation 

In this section two simulations of the complete model (9) are presented. The first simula- 
tion describes the behavior of the glucose-insulin regulatory system in the physiological 
case. The second simulation particularly shows the adaption of /3-cell mass under long 
term glucose infusion [12]. 

The physiological case - given in Figure 5 - has been simulated over 120 minutes with 
a high initial glucose value resulting for example from recent food intake. Parameters 
and initial values of the variables are given in Tables 1 and 2, respectively. The following 
description discusses the subplots. 

• Figure 5a: With the given parameters the threshold value for the cell cycle is 
G = 80 10 "^ nl - Glucose concentration G is above this level for 120 minutes. 
Therefore, the cell cycle reacts in the following way: Glucose values above the 
threshold increase the transition rate from G\- to S-phase. As the /3-cell cycle is a 
slow process, in the first 120 minutes an only slight increase in S- and G2/M-phase is 
detectable while the cell number in Gi-phase is decreasing. Increase in yS-cell mass, 
i.e., G\, takes more than 120 minutes. As glucose concentration is almost at the 
steady state G = 80 tggfgr towards the end of the simulation, the system will regulate 
itself without significant adaption of /3-cell mass. This is expected in the physiological 
case without abnormal exposure to glucose. 

• Figure 5b: The releasable amount of insulin X\ shows a biphasic behavior. There is a 
first peak release of stored insulin molecules and a second phase in consequence of 
provision of further insulin. 



x 10 










G 1 phase 




S phase 




G^/M phase 







time t in minutes 



time tin minutes 




time t in minutes 



time t in minutes 



Figure 5 Solution of the complete model over 120 minutes. The simulation presents the physiological 
case of the glucose-insulin regulatory system over 120 minutes after a high initial glucose value. The initial 
values for this simulation (with corresponding units) are given in Table 2 and the parameters in Table 1 . The 
solution of the complete model (9) was achieved numerically using Matlab ODE45. 
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• Figure 5c: With decreasing glucose concentration there are more packets with 
threshold value above the actual glucose level. For this reason and due to enhanced 
insulin provision P the amount of non-releasable insulin X2 increases. 

• Figure 5d: Glucose concentration G is decreasing from the high initial value as there 
is an increased concentration of insulin / in the blood. 

• Figure 5e: Blood insulin concentration I follows with some delay the releasable 
amount of insulin. It also shows the characteristic biphasic behavior of insulin release. 

• Figure 5f: Provision factor P shows an increase in presence of high glucose values 
and decreases as blood glucose decreases. 

The second simulation is done according to the experiments of Bonner-Weir et al. [12]. 
There, rats are given a high glucose infusion for 96 hours. After this time a significant 
increase in /3-cell mass is observable. The design of the experiment is assigned to the 
mathematical model (9) in the following way: the glucose production rate was increased 
from p 6 = 0.3 100 ^f min up to p 6 = 8.6806 100 mf m i n t0 account for a high concentrated 
glucose infusion. The resulting plots of the solution of model (9) are shown in Figure 6. 
Most important, a significant increase of /J-cells G\ can be seen in Figure 6a. This results 
from the persisting and severe hyperglycemia (Figure 6d) due to the high glucose produc- 
tion ps- With longer time of simulation the system will reach the stable steady state F^. 
The transfer of the experimental design in [12] to our model shows qualitatively that the 
glucose-insulin regulatory system is able to achieve euglycemia through adaption of /3-cell 
mass as it is stated in several publications (e.g. [6-8,13]). Note that the model currently 
assigns increase of /3-cell biomass exclusively to cell number (hyperplasia). It disregards 
an increase in cell size (hypertrophy) which occurs in [12]. This simplification which can 







G 1 phase 

S phase 

phase 











1000 2000 3000 4000 5000 

time t in minutes 




1000 2000 3000 4000 

time t in minutes 



1000 2000 3000 4000 

time t in minutes 




2000 3000 4000 

time t in minutes 



2000 3000 4000 

time t in minutes 



2000 3000 4000 

time t in minutes 



Figure 6 Solution of the complete model over 96 hours. The simulation presents an experimental 
situation with high glucose infusion over 96h.The parameter pe of glucose production was increased up to 
P6 = 8.6806 1Q0 ^ m|n . The fi-ceW mass increases due to the persisting hyperglycemia. With the adaption of 
fi- cell mass, seen in Figure 6a) the glucose-insulin regulatory system is able to reach euglycemia, i.e. the 
stable steady state F*. Other parameters and initial values for this simulation are given in Tables 1 and 2, 
respectively. The solution of the complete model (9) was achieved numerically using Matlab ODE45. 
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be overcome in further development of the model leads to an overestimation of changes 
in cell division rate as reflected in Figure 6a. Probably especially the first responses to 
elevated glucose concentrations are concerned. 

Discussion 

This work presents a mathematical model that describes three different negative regula- 
tion feedback loops of the glucose-insulin regulatory system: 

1. Immediate release of stored insulin molecules. 

2. Enhancement of provision of new insulin. 

3. Adaption of the fi-ceW cycle to metabolic demands. 

This is possible by incorporating an insulin secretion model describing storage and 
release of insulin molecules on the one hand and insulin provision on the other hand [4]. 
Furthermore, insulin provision is glucose dependent which allows its adaption to spe- 
cific demands of the organism at every time step. The third feedback loop is modeled via 
incorporation of the /3-cell cycle with glucose regulating the replication rate of the cells 
[11,13]. 

Several models of the glucose-insulin regulatory system, as e.g., [19-21], describe glu- 
cose, insulin and /J-cell mass dynamics, whereas our model shows the connection of 
glucose, and insulin concentrations with the /3-cell cycle as the main aspect. In this way 
the important role of glucose as regulator of the cell cycle [13] and the capability of the fi- 
cell mass to adapt to metabolic demands can be analyzed in detail. Hereby, the adaption 
of /J-cell mass is assigned exclusively to hyperplasia and disregards hypertrophy. 

The model conserves typical characteristics of the glucose-insulin regulatory system. 
The plots of the complete model in Figure 5 show biphasic insulin release represented 
through the biphasic shape of the releasable amount of insulin X\. This is a typical 
behavior of insulin release reported in several biological publications (e.g. [47,50]). 

Modeling insulin secretion based on [4] incorporates three feedback loops consisting of 
stored insulin, provision of further insulin, and variable /3-cell mass. Our model expands 
classic insulin secretion models (e.g. [4,44-46]) by a connection to the /J-cell cycle. 

The qualitative behavior of the model is illustrated with simulations. In the physiolog- 
ical case, shown in Figure 5, the /S-cell mass is sufficient to produce and release enough 
insulin to decrease glucose concentration and maintain euglycemia. With a second sim- 
ulation the adaption of the ^ -cell mass to increasing metabolic demands is presented. 
This situation occurs in long term studies with persisting hyperglycemia as it can be seen 
in Figure 6. In a simulation similar to the experiment in [12] increase of /J-cell mass via 
hyperplasia in 96 hours of hyperglycemia is shown. 

The analysis of the model gives the existence of two steady states. One describes a triv- 
ial fixed point with /J-cell numbers in all three phases equal to zero. The second steady 
state is a stable fixed point resulting from successfully achieving euglycemia. The trivial 
steady state is unstable while the stable non-trivial steady state will be attained in both sit- 
uations shown in Figures 5 and 6 with longer simulation times. The stability of the steady 
states is dependent on the underlying parameter values. As the parameters in our model 
are chosen in a way to describe a physiological situation the trivial steady state is unsta- 
ble and will not be reached. However, there are scenarios where the /J-cells in the model 
eventually die out. This could be due to an abnormally high apoptosis rate p± or artificially 
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holding of the glucose concentration below the threshold value G by external insulin infu- 
sion, for example. While a reduction of the replication rate due to hypoglycemia is shown 
in [11] a complete disappearance of the /3-cell mass is an implausible result. The death of 
the organism would happen prior to the extinction of the /3-cells. 

Our system allows for simulation of the glucose-insulin regulatory system assuming 
in vivo situations. Our model builds a theoretical basis for description and explana- 
tion of dynamics derived from biological experiments. It supports the understanding of 
metabolic processes. Biological assumptions can be verified and quantification of data 
and parameters can be achieved. Additionally, the model promotes understanding of the 
interplay of the three different regulation feedback loops. The model is able to describe 
metabolic dynamics of the glucose-insulin regulatory system also for the pathological case 
of type 1 or type 2 diabetes. 

To illustrate one possible modification of the system a type 2 diabetes-like simulation is 
done. Type 2 diabetes is characterized by insulin resistance of target cells, mainly muscle 
and fat cells. In consequence, these cells are not able to take up enough glucose from the 
blood. In this case the insulin sensitivity of the body cells is down-regulated. To simulate 
this situation the model parameter for insulin sensitivity pg is decreased arbitrary from 
the value pg = 360 x 10~ 3 to ps = 360 x 10~ 5 while the other parameters given in Table 1 
stay the same. This modification corresponds to a lower reaction of the target cells to 
insulin and therefore a decreased uptake of glucose from the blood. 

Figure 7 shows that the blood glucose concentration G in the pathological case of insulin 
resistance (dashed line) decreases slower than in the physiological case (solid line). The 
body cells take up less glucose from the blood and therefore the hyperglycemia lasts 
longer in the pathological case. 

There are also possibilities to simulate the regulatory system in a type 1 diabetes sce- 
nario. It can be done for example by increasing the apoptosis rate pn in the cell cycle 
model. This results in dying /3-cells which is characteristic for this autoimmune disease. 
For a more detailed discussion about /i-cell mass in a type 1 diabetes scenario see Klinke 
[51]. There, /S-cell mass at onset of type 1 diabetes is concerned depending on body weight 
and the patient's age. 
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Figure 7 Insulin resistance of target cells. Blood glucose concentration at different values of insulin 
sensitivity pg. Physiological value pg = 360 x 1 0~ 3 (solid line) and pathological value pg = 360 x 1 0~ 5 
(dashed line). 
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In summary, our model is a basic approach to understand the processes within the 
glucose-insulin regulatory system connected to the /i-cell cycle. It offers a wide range 
of possible modifications to incorporate further processes and can be adapted to many 
biological questions. 

Appendix 

Derivation of glucose dependent insulin storage dynamics 

With time dependent glucose the major difference to the original model of Grodsky 
consists in time dependence of upper and lower integral bounds for X\ and Xi: 

d I f G(t) 
X 1 (t)=-[l W,L)d»], 



f 

Jo 

= "ff 



x ^) = —a / mt)de\. 

>G(t) / 

To determine these integrals, the differential equation for f (6, t), given in [4], - has to be 
solved first. Then this solution can be applied to the expressions for X\ and Xi. In a final 
step we differentiate these terms with respect to time t to achieve differential equations 
X\ and^2. 

Solution of the differential equation £(0,f) 

As only the packets with threshold 9 below glucose concentration G are relevant for 
insulin secretion, two cases have to be distinguished. There is a differential equation for 

0 < G(t), 

1(0, t) = -p9 t) +fy(0)P(t) -fX max m t) +ho(B)Xoo(t), (17) 
and for 6 > G(t), 

k(B, t) =f&(6)P(t) -fX max m t) +ho(B)Xoo(t), (18) 

with 

poo POO 

MO) = HO,0),Xoo(t)= / W'.W =Xi(t)+X 2 (t)andX max = / $ 0 (6')d6'. 

Jo Jo 

The only difference between the two equations is the term for insulin secretion, 
— p9 |(0, t), that occurs in the first but not in the second equation. Therefore, we restrict 
the following exploration to the case 9 < G(t). 

The system that has to be solved is a linear non-homogeneous ordinary differential 
equation for which there are standard solution methods like variation of constants. Using 
this method, the set of solutions for differential equation (17) is given as 



f(6U) = j e -(P9+/*«,*)f jf e (P9+fx max )T j^p (T) + ^ Xoo(T) j £ o(6 i)j T + C 
and for the differential equation (18) as 



e -fXm 



t jf /w[jP(r)+/I M (r)]? 0 (Wr + C 



Ce 



The constant C incorporates the initial conditions by C = %o(9). 
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Derivation of the differential equations Xi (f) andX2(f) 

To derive a differential equation for the total amount of releasable insulin at glucose 
concentration G the expression 

pG(t) 

Xi(t) = / t-{0,t)dO 
Jo 

has to be differentiated with respect to time t. Attaching the solution of %(0,t) to X\ 
results in the following expression: 

pG(t) r ft _ 

Xi(f) = e- (p9+fXmax)t J $o(6)d6 J e (P9+ - fXmax)r \fP{r) +/X oc (t)J dr + 1 
The differentiation of this equation with respect to time t then results in 
^-Xit) = * lit) ^(G(ihi)G(i) -p 9 X x {t) 

fG(t) _ poo 

+// s*(o,t)doix 1 (t)+x 2 (t)]-f / $*{e,t)dex x (t) 

Jo Jo 



Jo 



G(t) 

t(O,t)d0P(t), 



with target distribution f * (0, t). 

Analogously, we derive the differential equation for the total amount of non-releasable 
insulin at glucose concentration G 

poo 

x 2 (t) = / H0,t)d6. 

JG(t) 

The differential equation is given as 

£x 2 (f)=X 2 (t) = ~™\ J *(G(t),t)G(t) 
dt j m %*(6,t)d6 

/•OO pOO 

+// $*(0,t)d0[X 1 (t)+X 2 (t)]-f ?{0,t)dOX 2 (t) 

JG(t) Jo 
poo 

+// t(O,t)d0P(t), 

JG(t) 

with target distribution t). 
Transition functions 

The transition functions u\,...,u\ of the compartment model for insulin secretion are 
complex expressions with integrals over the target distribution in the insulin storage. With 
the concrete packet distribution given in [4] the integrals can be determined explicitly in 
terms of Hill functions: 

Git) 



. m , Wfl X max (t) G{t) k Gi(t) G(t) k 

? (9, t)d9 = —7 ——7 = X„ 



/ s " C k + Guv ' ""'"' v G'i(0) e ■!- GUV ' 

0 



and 



X max (t) C k - Gi(t) C k 
f * (9, t)d9 = , , = X„ 



J b ' C k \ G(L) k ' ' M " nA 'Gi(0) C<- +GU) k ' 

G(t) 



The parameters / and / are proportionality factors. The constant character of these fac- 
tors has to be changed in our case due to variable glucose and the connection to the /J-cell 
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cycle that influences the maximum amount of insulin per pancreas. Therefore, / is chosen 
glucose dependent and with a condition to ensure that 

X l (t)+X 2 (t)<X max (t). (19) 
The function / has the form 

p 9 f C k G(t) h ~ k 



/(*) = 



Xmax§$ [fG(t) h (l + C k G(tT k ) - P9 (P h Q + G(t) h )] ' 



This expression was found by analyzing condition (19) in steady state situation. 
With these preliminaries the transition functions can be given explicitly: 

-p 9 f G(t)*~* C 2k 



u\(t) 



u 2 (t) = 



[fG(t)Hl + C k G(tr k ) - P9 (P h 0 + G(t) h )] [C k + G(t) k ] ' 

-p 9 f G(t) h C k 

[fG(t) h (l + C k G(t)~ k ) - P9 {P h Q + G(t)*)] [C* + G(t) k ] ' 
1 G(t) k 



" 3(t) -/Jfm8 *Gi(0) C +GU)^ 
114(f) =fX max - 



1 



Gi(0) C* + G(t)*' 
These expressions were used for the simulations in Figures 5 and 6. 
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